function Me = findMachNumber(pratio,gamma)
% function with input of pressure ratio and gamma
% output mach number
% finds mach number for a given pressure ratio by solving implicit equation
% as p/pa = (1+(gamma-1)/2*Me^2)^(gamma*(gamma-1))
% this is the isentropic flow mach number at a converging throat

data.pratio = pratio;
data.gamma = gamma;

% Me = fzero(@machEqn,0.5,optimset(),data);
Me = newtonZero(@machEqn,.5,data,1e-5);


end


function [residual,slope] = machEqn(M,data)

residual = data.pratio - (1+(data.gamma - 1)/2*(M^2))^...
    (data.gamma*(data.gamma-1));

slope = -data.gamma*(data.gamma-1)*...
    ((1+(data.gamma-1)/2*(M^2))^(data.gamma*(data.gamma-1)-1))*...
    ((data.gamma-1)*M);
end


function x = newtonZero(fhandle,x0,data,tol)

[r,m] = feval(fhandle,x0,data);
x = x0;

stabFact = 0.5;
maxIts = 10000;

itN = 0;
while abs(r) > tol
    itN = itN + 1;
    if itN > maxIts
        error('exceeded maxits in newton')
    end
    x0 = x;
    r0 = r;
    m0 = m;
    
    x = x0 - stabFact*r0/m0;
    
    [r,m] = feval(fhandle,x0,data);
end
end
